
/* THIS DO-FILE WAS WRITTEN IN MAY 2019 BY MARGHERITA COMOLA,
IT RUNS ALL RESULTS FOR THE EMPIRICAL ILLUSTRATION ON NEPALESE DATA
(`Treatment Effects Accounting for Network Changes', Comola Prina, REStat) */

clear all
set maxvar 30000
set matsize 11000
set more off

/*******************************************************************
********************************************************************
*** COMPUTE AND SAVE THE INTERACTION MATRICES
*******************************************************************
********************************************************************/

use dyadic_dataset_ComolaPrina.dta, clear

sort HHNO1 HHNO2
fillin  HHNO1 HHNO2
replace link_bin0=0 if link_bin0==.
replace link_bin1=0 if link_bin1==.

mata 

/* binary matrix */

G0 = st_data(., "link_bin0")
G0 = colshape(G0, 915) 

G1 = st_data(., "link_bin1")
G1 = colshape(G1, 915) 

/* semi row-standardization */

G0_S=rowsum(G0,.)
G0=G0:/G0_S

G1_S=rowsum(G1,.)
G1=G1:/G1_S

mata drop G0_S G1_S

/* missing values */

G0=editmissing(G0, 0)
G1=editmissing(G1, 0)

/* delta_G */

deltaG=G1-G0

st_matrix("G0", G0)
st_matrix("G1", G1)

end

/* re-import into Stata */

clear
svmat G0
stack G01-G0915, into (link_s0)
gen serial=_n
sort serial 
drop _stack
save temp0.dta, replace 

clear
svmat G1
stack G11-G1915, into (link_s1)
gen serial=_n
sort serial 
drop _stack
save temp1.dta, replace 

/* merge into the dyadic data */

use dyadic_dataset_ComolaPrina.dta, clear

fillin HHNO1 HHNO2
sort  HHNO1 HHNO2
gen serial=_n

sort serial
merge serial using temp0
drop _merge
capture erase temp0.dta

sort serial
merge serial using temp1
drop _merge
capture erase temp1.dta

keep if _fillin==0
drop _fillin serial
replace link_s0=. if HHNO1==HHNO2
replace link_s1=. if HHNO1==HHNO2

/*******************************************************************
********************************************************************
*** DYADIC ANALYSIS
*******************************************************************
********************************************************************/

/* des stats*/

gen some_treated=(ITT_H11==1 | ITT_H21==1)
gen both_treated=(ITT_H11==1 & ITT_H21==1)
gen one_treated=((ITT_H11==1 & ITT_H21==0) | (ITT_H11==0 & ITT_H21==1))

tabstat link_bin0 link_s0 both_married0 abs_diff_chl160 abs_diff_HHmembers0 some_shock_livestock0 some_shock_death0 if HHNO1!=HHNO2, stats (count mean min max sd) columns(statistics) longstub varwidth(30) labelwidth(30)
tabstat link_bin1 link_s1 both_married1 abs_diff_chl161 abs_diff_HHmembers1 some_shock_livestock1 some_shock_death1 some_treated both_treated one_treated if HHNO1!=HHNO2, stats (count mean min max sd) columns(statistics) longstub varwidth(30) labelwidth(30)

tab link_bin0 link_bin1, row
tab link_bin0 link_bin1 if some_treated==0, row 
tab link_bin0 link_bin1 if some_treated==1, row

/* format */

foreach var in link_s link_bin abs_diff_HHmembers abs_diff_chl16 both_married some_shock_livestock some_shock_death { 
gen delta_`var'= `var'1-`var'0
drop `var'1 `var'0
}

/* dyadic regressions */

global controls = "delta_both_married delta_abs_diff_chl16 delta_abs_diff_HHmembers delta_some_shock_livestock delta_some_shock_death"

regress delta_link_s some_treated, noc 
regress delta_link_s some_treated $controls

matrix coeffs_dyadic=e(b)
matrix der=coeffs_dyadic[1,1]

mata: der=st_matrix("der")

/* `dyadic' command: dyadic s.e. adjusted for incomplete network data */
   
dyadic delta_link_s some_treated, noc id(HHNO1 HHNO2) inc
dyadic delta_link_s some_treated $controls, id(HHNO1 HHNO2) inc

/*******************************************************************
********************************************************************
 PEER EFFECT REGRESSIONS:
 - the dependent variable is meat consumption
 - semi row-standardized interaction matrix 
 - 2nd and 3rd order IVs 
 - clustered boostrapped s.e.
*******************************************************************
*******************************************************************/

*******************************************************************
* COMPUTE THE RELEVANT MATRICES
*******************************************************************

use individual_dataset_ComolaPrina.dta, clear

mata

/* save individual characteristics */

ITT= st_data(., "ITT")
y0= st_data(., "y0")
y1= st_data(., "y1")

/* vectors of interest */

G0_ITT=G0*ITT
G0_y0=G0*y0
G0_y1=G0*y1
G0_deltay=G0*(y1-y0)
deltaG_ITT=deltaG*ITT
deltaG_y1=deltaG*y1

/* IVs of second order */

IV_1=G0*G0*ITT
IV_2=deltaG*deltaG*ITT
IV_3=G0*deltaG*ITT
IV_4=deltaG*G0*ITT

/* IVs of third order */

IV_5=G0*G0*G0*ITT
IV_6=deltaG*deltaG*deltaG*ITT
IV_7=G0*G0*deltaG*ITT
IV_8=deltaG*deltaG*G0*ITT
IV_9=G0*deltaG*deltaG*ITT
IV_10=deltaG*G0*G0*ITT
IV_11=G0*deltaG*G0*ITT
IV_12=deltaG*G0*deltaG*ITT

results=(G0_ITT, G0_y0, G0_y1, G0_deltay, deltaG_ITT, deltaG_y1, IV_1, IV_2, IV_3, IV_4, IV_5, IV_6, IV_7, IV_8, IV_9, IV_10, IV_11, IV_12)

st_matrix("result", results)

end

/*re-import into Stata*/

svmat result

ren result1 G0_ITT
ren result2 G0_y0
ren result3 G0_y1
ren result4 G0_deltay
ren result5 deltaG_ITT
ren result6 deltaG_y1
ren result7 IV_1
ren result8 IV_2
ren result9 IV_3 
ren result10 IV_4 
ren result11 IV_5 
ren result12 IV_6 
ren result13 IV_7 
ren result14 IV_8 
ren result15 IV_9 
ren result16 IV_10 
ren result17 IV_11
ren result18 IV_12

/* village size */

preserve

bysort village_code: gen village_tot=_N
bysort village_code: gen serial=_n
drop if serial>1
sort HHNO
keep village_tot

mata
SIZE= st_data(., "village_tot")
end 

restore

*******************************************************************************
* PEER EFFECT REGRESSIONS 
*******************************************************************************

cap erase results_pe.txt
cap erase results_pe.xml

set seed 07031979

/* no PE*/ 

regress delta_y ITT, nocons vce (boot, r(100) cluster(village_code))
outreg2 using results_pe, excel replace dec(2) ctitle(no PE)

/* standard PE, IV */

ivregress 2sls delta_y ITT (G0_deltay = IV_1-IV_12) G0_ITT, nocons vce (boot, r(100) cluster(village_code))
outreg2 using results_pe, excel append dec(2) ctitle("st PE, IV")

matrix coeffs_stat=e(b)
matrix beta1_stat= coeffs_stat[1,1]
matrix gamma_stat= coeffs_stat[1,2]
matrix delta1_stat= coeffs_stat[1,3]

mata: beta1_stat=st_matrix("beta1_stat")
mata: gamma_stat=st_matrix("gamma_stat")
mata: delta1_stat=st_matrix("delta1_stat")

/* dynamic PE, IV*/

ivregress 2sls delta_y ITT (G0_deltay deltaG_y1 = IV_1-IV_12) G0_ITT deltaG_ITT, nocons vce (boot, r(100) cluster(village_code))
outreg2 using results_pe, excel append dec(2) ctitle("dyn PE, IV")

matrix coeffs_dyn=e(b)
matrix  beta1_dyn= coeffs_dyn[1,1]
matrix  beta2_dyn= coeffs_dyn[1,2]
matrix gamma_dyn= coeffs_dyn[1,3]
matrix  delta1_dyn= coeffs_dyn[1,4]
matrix  delta2_dyn= coeffs_dyn[1,5]

mata: beta1_dyn=st_matrix("beta1_dyn")
mata: beta2_dyn=st_matrix("beta2_dyn")
mata: gamma_dyn=st_matrix("gamma_dyn")
mata: delta1_dyn=st_matrix("delta1_dyn")
mata: delta2_dyn=st_matrix("delta2_dyn")

test G0_deltay=deltaG_y1
cap erase results_pe.txt

********************************************************************
********************************************************************
*** DISTRIBUTION OF INDIVIDUAL EFFECTS
********************************************************************
********************************************************************

/* individual panel */

preserve

keep HHNO village_code y0 G0_y0

ren y0 yh
ren G0_y0 G0h_y

gen deltaGh_y=0
gen ITTh=0
gen G0h_ITTh=0
gen deltaGh_ITTh=0

forvalues x=1/12 {
gen IV_`x'=0
}

gen time=0
save temp1.dta, replace 

restore 

preserve

ren y1 yh
ren G0_y1 G0h_y
ren deltaG_y1 deltaGh_y
ren ITT ITTh
ren G0_ITT G0h_ITTh
ren deltaG_ITT deltaGh_ITTh

keep HHNO village_code yh G0h_y deltaGh_y ITTh G0h_ITTh deltaGh_ITTh IV_1-IV_12
gen time=1

append using temp1.dta
capture erase temp1.dta

xtset HHNO time

/* save the individual effects */

qui tab HHNO, gen (code)
xi: qui ivregress 2sls yh ITTh (G0h_y deltaGh_y = IV_1-IV_12) G0h_ITTh deltaGh_ITTh code2-code915, nocons
matrix miu=e(b)
matrix miu= miu[1,5...]'

clear
svmat miu, name(miu)
replace miu1=0 in 1

gen serial=_n
sort serial
save miu.dta, replace 

restore

gen serial=_n
sort serial
merge serial using miu.dta
drop _merge serial
capture erase miu.dta

********************************************************************
********************************************************************
*** TOTAL TREATMENT EFFECT
********************************************************************
********************************************************************

/******************************
* STATIC MODEL 
******************************/

mata 

/* full matrix of cross derivatives */ 

A1 = J(915, 915, . )
for (i=1; i<=915; i++) {
e= J(915,1,0)
e[i]=1
A1[.,i]=(luinv(I(915)-beta1_stat*G0))*(gamma_stat*e+delta1_stat*G0*e)
}

/* direct effect */
 
D1 = diagonal(A1)
DIRECT_static=mean(D1)
DIRECT_static

/* total effect */ 

C1 = colsum(A1)'
TOTAL_static=mean(C1)
TOTAL_static

/* indirect effect */ 

S1=(C1-D1)
INDIRECT_static=mean(S1)
INDIRECT_static

end

/************************
 DYNAMIC MODEL, NO MIUs
************************/

mata 

V1=J(60,60,der)
V2=J(26,26,der)
V3=J(74,74,der)
V4=J(82,82,der)
V5=J(28,28,der)
V6=J(119,119,der)
V7=J(26,26,der)
V8=J(47,47,der)
V9=J(38,38,der)
V10=J(48,48,der)
V11=J(25,25,der)
V12=J(12,12,der)
V13=J(51,51,der)
V14=J(36,36,der)
V15=J(74,74,der)
V16=J(33,33,der)
V17=J(61,61,der)
V18=J(11,11,der)
V19=J(64,64,der)

BLOCK=blockdiag(V1,blockdiag(V2,blockdiag(V3,blockdiag(V4,blockdiag(V5,blockdiag(V6,blockdiag(V7,blockdiag(V8,blockdiag(V9,blockdiag(V10,blockdiag(V11,blockdiag(V12,blockdiag(V13,blockdiag(V14,blockdiag(V15,blockdiag(V16,blockdiag(V17,blockdiag(V18,V19))))))))))))))))))

/* full matrix of cross derivatives */

A2 = J(915, 915, . )

for (i=1; i<=915; i++) {

e= J(915,1,0)
e[i]=1

/* derivative of G wrt ITT_k */

DER_G=J(915, 915, 0)
DER_G[i,.]=BLOCK[i,.]
DER_G[.,i]=BLOCK[.,i]

S=luinv(I(915)-beta1_dyn*G0-beta2_dyn*deltaG)

M=gamma_dyn*ITT+delta1_dyn*G0*ITT+delta2_dyn*deltaG*ITT

DER_S=S*beta2_dyn*DER_G*S

DER_M=gamma_dyn*e+delta1_dyn*G0*e+delta2_dyn*DER_G*ITT+delta2_dyn*deltaG*e

A2[.,i]=DER_S*M+S*DER_M

}

/* direct effect */ 

D2 = diagonal(A2)
DIRECT_dynamic=mean(D2)
DIRECT_dynamic

/* total effect */ 

C2 = colsum(A2)'
TOTAL_dynamic=mean(C2)
TOTAL_dynamic

/* indirect effect */ 

S2=(C2-D2)
INDIRECT_dynamic=mean(S2)
INDIRECT_dynamic

end
 
/*****************************
 DYNAMIC MODEL, ESTIMATED MIUs
*****************************/

mata 

miu= st_data(., " miu1")

V1=J(60,60,der)
V2=J(26,26,der)
V3=J(74,74,der)
V4=J(82,82,der)
V5=J(28,28,der)
V6=J(119,119,der)
V7=J(26,26,der)
V8=J(47,47,der)
V9=J(38,38,der)
V10=J(48,48,der)
V11=J(25,25,der)
V12=J(12,12,der)
V13=J(51,51,der)
V14=J(36,36,der)
V15=J(74,74,der)
V16=J(33,33,der)
V17=J(61,61,der)
V18=J(11,11,der)
V19=J(64,64,der)

BLOCK=blockdiag(V1,blockdiag(V2,blockdiag(V3,blockdiag(V4,blockdiag(V5,blockdiag(V6,blockdiag(V7,blockdiag(V8,blockdiag(V9,blockdiag(V10,blockdiag(V11,blockdiag(V12,blockdiag(V13,blockdiag(V14,blockdiag(V15,blockdiag(V16,blockdiag(V17,blockdiag(V18,V19))))))))))))))))))

/* full matrix of cross derivatives */

A2 = J(915, 915, . )

for (i=1; i<=915; i++) {

e= J(915,1,0)
e[i]=1

DER_G=J(915, 915, 0)
DER_G[i,.]=BLOCK[i,.]
DER_G[.,i]=BLOCK[.,i]

S=luinv(I(915)-beta1_dyn*G0-beta2_dyn*deltaG)

M=gamma_dyn*ITT+delta1_dyn*G0*ITT+delta2_dyn*deltaG*ITT+miu

DER_S=S*beta2_dyn*DER_G*S

DER_M=gamma_dyn*e+delta1_dyn*G0*e+delta2_dyn*DER_G*ITT+delta2_dyn*deltaG*e

A2[.,i]=DER_S*M+S*DER_M

}

/* direct effect */ 

D2 = diagonal(A2)
DIRECT_dynamic=mean(D2)
DIRECT_dynamic

/* total effect */ 

C2 = colsum(A2)'
TOTAL_dynamic=mean(C2)
TOTAL_dynamic

/* indirect effect */ 

S2=(C2-D2)
INDIRECT_dynamic=mean(S2)
INDIRECT_dynamic

end

/*******************************************************************************
CHECK FOR WEAK IDENTIFICATION (Beugnot et al. EER 2019) 
*******************************************************************************/

mata

M1 = I(915)
M1 = colshape(M1, 1) 

M2 = G0
M2 = colshape(M2, 1) 

M3 = deltaG
M3 = colshape(M3, 1) 

M4 = G0*G0
M4 = colshape(M4, 1) 

M5 = deltaG*deltaG
M5 = colshape(M5, 1) 

M6 = G0*deltaG
M6 = colshape(M6, 1) 

M7 = deltaG*G0
M7 = colshape(M7, 1) 

condition_matrix=(M1, M2, M3, M4, M5, M6, M7)
condition_number=cond(condition_matrix)

condition_number

end

